The URDME manual 

Version 1.0 



Josef CuUhed Stefan Engblom' Andreas Hellander^ 

December 15, 2008 

^Div of Scientific Computing, Dept of Information Technology 
Uppsala University, P. 0. Box 337, SE-75105 Uppsala, Sweden 
email: stefane@it.uu.se, andreas.hellaiider@it.uu.se 

Abstract 

We have developed URDME, a general software for simulation of 
stochastic reaction-diffusion processes on unstructured meshes. This 
allows for a more flexible handling of complicated geometries and 
curved boundaries compared to simulations on structured, Cartesian 
meshes. The underlying algorithm is the next subvolume method 
(NSM), extended to unstructured meshes by obtaining diffusion jump 
coefficients from the finite element formulation of the corresponding 
macroscopic equation. 

In this manual, we describe how to use the software together with 
Comsol Multiphysics 3.4 and Matlab to set up simulations. We 
provide a detailed account of the code structure and of the available 
interfaces. This makes modifications and extensions of the code pos- 
sible. We also give two detailed examples, in which we describe the 
process of simulating and visualizing two models from the systems 
biology literature in a step-by-step manner. 

Keywords: URDME, reaction-diffusion master equation, stochastic 
chemical kinetics, stochastic spatial simulation, unstructured meshes. 



1 Introduction 

Stochastic simulation methods are frequently used to study the behavior of 
cellular control systems modeled as continuous-time discrete-space Markov 
processes (CTMC). Compared to the most frequently used deterministic 



1 



model, the reaction rate equations, the mesoscopic stochastic description 
can capture effects from intrinsic noise on the behavior of the networks 

[1, 6, 23, 24, 27]. 

In the discrete mesoscopic model the state of the system is the copy 
number of the different chemical species and the reactions are usually as- 
sumed to take place in a well-stirred reaction volume. The chemical master 
equation is the governing equation for the probability density, and for small 
to medium sized systems it can be solved by direct, deterministic methods 
[8, 9, 13, 19, 22, 26]. For larger models however, exact or approximate kinetic 
Monte Carlo methods [15, 16] are frequently used to generate realizations of 
the stochastic process. Many different hybrid and multiscale methods have 
also emerged that deal with the typical stiffness of biochemical reactions 
networks in different ways, see e.g. [2, 4, 17, 20, 21, 25]. 

Many processes inside the living cell can not be expected to be explained 
in a well-stirred context. The natural macroscopic model is the reaction- 
diffusion equation and it has the same limitations as the reaction rate equa- 
tions. By discretizing the space with small subvolumes it is possible to model 
the reaction-diffusion process by a CTMC in the same formalism as for the 
well-stirred case. A diffusion event is now modeled as a first order reaction 
from a subvolume to an adjacent one and the state of the system is the 
number of molecules of each species in each subvolume. The corresponding 
master equation is called the reaction-diffusion master equation (RDME) 
and due to the very high dimensionality it cannot be solved by deterministic 
methods for realistic problem sizes. 

The RDME has been used to study biochemical systems in [5, 12]. Here 
the next subvolume method (NSM) [5] , an extension of Gibson and Bruck's 
next reaction method (NRM) [14], was suggested as an efficient method for 
realizing sample trajectories. An implementation on a structured Cartesian 
grid is freely available in the software MesoRD [18]. 

The method was extended to unstructured meshes in [10] by making 
connections to the finite element method (FEM). This has several advantages, 
but the most notable one is the ability to handle complicated geometries in a 
flexible way which is particularly important when internal structures of the 
cell must be taken into account. 

This manual describes the software URDME which implements the un- 
structured extension of NSM as suggested in [10] . The purpose with the code 
is to provide an efficient, modular implementation that is easy to use for sim- 
ulating and studying a particular model in an applied context, but also for 
developing and testing new approximate methods. We achieve this by relying 
on commercial software for the geometry definition, meshing, preprocessing 
and visualization and use a highly efficient computational core written in 
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Ansi C. This keeps the implementation of the actual Monte Carlo simula- 
tion small and easily extendible, while the user benefits from the advanced 
pre- and postprocessing capabilities of modern FEM software. In URDME, 
we have chosen to provide an interface to Comsol Multiphysics 3.4 [3]. 

The rest of this manual is organized as follows. In Section 2 we recapit- 
ulate the mesoscopic reaction-diffusion model and show how the stochastic 
diffusion intensities are obtained from a FEM discretization of the diffusion 
equation. An overview of the code structure is also offered in Section 3. The 
details concerning the input to the code, the provided interface to Comsol 
Multiphysics 3.4 and the way models should be specified are found in 
Section 4. Finally, two models are set up and simulated in a step- by-step 
manner in Section 5, which is followed in Section 6 by a short discussion of 
performance. 

2 Background 

In this section we briefiy describe how reaction and diffusion events are mod- 
eled and how we obtain the diffusion rate constants when the domain is 
discretized using an unstructured mesh. For a fuller introduction to the 
subject along with many additional references, consult [7]. 

The computational core of URDME is based on the next subvolume 
method (NSM) [5], but it has been adapted to simulation on unstructured 
meshes by supporting a more general input format. Details concerning the 
actual simulation algorithms can be found in Appendix A. 

2.1 Mesoscopic chemical kinetics 

In a well-stirred chemical environment reactions are understood as transitions 
between the states of the integer- valued state space counting the number of 
molecules of each of D different species. The intensity of a transition is 
described by a reaction propensity defining the transition probability per 
unit of time for moving from the state x io x — N^-^ 

x^x-Nr, (2.1) 

where A^^ G Z-^ is the transition step and is the rth column in the stoichio- 
metric matrix N. Eq. (2.1) defines a continuous-time Markov chain over the 
positive D-dimensional integer lattice. 

When the reactions take place in a container of volume Q, it is sometimes 
useful to know that the propensities often satisfy the simple scaling law 

LOrix) = nUr{x/n) (2.2) 
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for some function Ur which does not involve Q,. Intensities of this form are 
called density dependent and arise naturally in a variety of situations [11, 
Ch. 11]. 

2.2 Mesoscopic diffusion 

In the mesoscale model, a diffusion event is modeled as a first order reaction 
taking species Si in subvolume Q from its present subvolume to an adjacent 
subvolume 

Su ^ (2.3) 

where x/j is the number of molecules of species / in subvolume i. On a uniform 
Cartesian mesh such as those used in MesoRD [18], the rate constant takes 
the value Uij — ^/h? where h is the side length of the subvolumes and 7 
is the diffusion constant. In URDME we use an unstructured mesh made 
up of tetrahedra and the rate constants are taken such that the expected 
value of the number of molecules divided by the volume (the concentration) 
converges to the solution obtained from a consistent FEM discretization of 
the diffusion equation 

ut = 7Aii. (2.4) 

Using piecewise linear Lagrange elements and mass lumping, we obtain the 
discrete problem 

Ut = M-^Ku (2.5) 

where M is the lumped mass matrix and K = {kij} is the stiffness matrix. 
The rate constants on the unstructured mesh are then given by 

ttij = 'o^^iji (2-6) 

where f2j is the diagonal entry of M and can be interpreted as the volume 
of the dual element associated with mesh node i (see Figure 2.1). For more 
details, consult [10]. 

The assumption made in the mesoscopic model is that molecules are well- 
stirred within a dual cell. These dual cells correspond to the cubes of the 
staggered grid in a Cartesian mesh. 

3 Code overview 

In this section we give an overview of the structure of the code. The compu- 
tational core of URDME is an efficient implementation of NSM and addition- 
ally consists of a small set of routines written in Ansi C, Matlab and Comsol 
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Figure 2.1: A 2D example of an unstructured triangular mesh. The primal 
mesh is shown in dashed and the dual in solid. Within each dual element 
the system is assumed to be well-stirred, and molecules can jump from each 
dual cell to the neighboring ones. 

Script. Table 3.1 shows the directory structure of URDME together with a 
short description of each routine. 

The software logically consists of three major parts. From bottom to top 
those are the Ansi C kernel, a Matlab/Mex-interface to the kernel and an 
interface to the FEM software. By design, the computational core is stand- 
alone, and in principle, for each FEM software a separate interface can be 
provided. This interface needs to export the internal representation of the 
problem and the mesh and convert it to the input format of rdme_solve 
which is detailed in Section 4. The actual call to rdme_solve is made via 
a gateway routine which could be a C mainO function or as in our case, 
through a Matlab/Mex construction. 

We have provided an interface to Comsol Multiphysics 3.4, which is 
used to specify the problem, create the mesh and provides for postprocessing. 
Since Comsol Script is very similar to the Matlab language, it is convenient 
to supply the gateway routine as a Mex-function using the Matlab C API. The 
Comsol Script routine f em2rdme assembles the matrix with diffusion rate 
constants D and the lumped mass matrix M. They are written to a . mat file 
together with information of the mesh and the different subdomains. This 
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Directory 


File(s) 


Description 


src 


rdme.c 
hinhpan c 
report.c 
mexrdme.c 


The computational core rdme_solve. 
The liinarv hean used in rdme solvs 
Report function used in rdme_solve. 
Mex-interface to rdme_solve. 


include 


rdme.h 

propensities. h. 
binheap.h 

i C^iJi t .11 


Header for rdme.c. 

Definition of the propensity function 
datatype. 

Header for binheap.c. 

I-Tpfirip'r TCw vt^T\rwi' o 

llCctLlCl liJl IC^VJlt.*^. 


comsol 


fem2rdme.m 
rdme2fem.m 


Comsol-script converting Comsol's FEM- 
struct to a valid rdme_solve-input. 
Comsol-script for conversion of the out- 
put of rdme.m to the solution format in 

ririTncinl Tnp "nnmo'^p to nnfain a A^alin 

V^WlilO W _L . i lie IJ HI IJWOC lo Uw yjVJhdLLL Ci VcHH.1 

FEM-struct so that e.g. Comsol's postpro- 
cessing can be used. 


msrc 


rdme.m 
make.m 
startup, m 


Mat lab- wrapper for mexrdme. Handles 

error checking of the input. 

Makefile for building the Mex-file 

mexrdme. 

Initiahzes URDME. 


examples 


(various) 


Contains files specifying the two examples 
studied in detail in Section 5. 


doc 


manual.pdf 


The most recent version of this manual. 



Table 3.1: Overview of the files and routines that make up URDME. 



.mat file can subsequently be loaded into Mat lab where some of the other 
inputs are specified and the simulation is made via rdme. 

The chemical reactions in the model can be specified in two different ways. 
Either as inline propensities specified from within Matlab or as compiled 
propensities where the rate functions are defined in a separate model file 
written in Ansi C. Inline propensities can be used to define basic polynomial 
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rate laws, while the compiled propensities are completely general. The details 
of this are found in Section 4. 

When the model has been completely specified, the call to rdme_solve is 
made via the interface pair rdme.m/mexrdme.c. 

At first, this construction may appear both cumbersome and inconve- 
nient. However, the structure is transparent and makes it possible to separate 
the computational core from the actual FEM software used and development 
and maintenance of the core is made easier. Also, there are inputs to the 
stochastic simulation that are hard to standardize such as initial conditions 
or "exotic" diffusion over manifolds. Using Matlab as an interface facilitates 
handling of input and non-standard pre- and postprocessing since both the 
input and output data are available in Matlab. Hence the user can provide 
postprocessing routines which is particularly valuable in experimental stages 
of model development. 

It is of course possible to eliminate the dependence on Matlab by pro- 
viding another gateway in the form of for example a C mainO routine, that 
directly reads the .mat file provided by the Comsol interface. This might be 
included in a later release. 

4 Details and specifications 

In this section we give a detailed description of the input to rdme_solve and 
we also give some details concerning the provided interface routines. Any 
gateway routine (see for example 'rdme.m' and 'mexrdme.c') should carefully 
ensure that the input has the correct format since the computational kernel 
puts very little effort in error checking. 

4.1 Input to rdiiie_solve 

Tables 4.1 and 4.2 summarize the input to rdme_solve. For precise type 
definitions, consult the preamble of the source file 'rdme.c'. 

The diffusion matrix D and the (diagonal of) the lumped mass matrix 
vol as well as the vector giving the (generalized) subdomain number of each 
subvolume, sd, are generated by the Comsol interface routine fem2rdme. 
How sd can be specified in the Comsol GUI is explained in Section 5. The 
generalized data vector data may be used to pass any additional arguments 
to the compiled propensity functions. 

The stoichiometry matrix N, the dependency graph G, the initial condition 
uO, the vector of output times tspan as well as the definition of the inline 
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propensities (if they are used) are explicitly specified in Matlab and passed 
on to rdme_solve via the Mex-interface in rdme, cf. Section 4.2.1. 

If compiled propensities are used, they are specified as a vector of func- 
tion pointers, prop. This is slightly more involved and is explained in Sec- 
tion 4.2.2. 

The above mentioned data is passed to sd via the Matlab interface rdme. 
The other inputs (see 'rdme.c') are implicitly defined by the data structures, 
and are obtained automatically in mexrdme, where also the report function 
report is defined. The functionality of the report function is stated in the 
help section of 'rdme.m' and can easily be modified to suit the need of the 
user. 

4.2 Specifying propensities for chemical reactions 

We have provided two separate methods to specify the reaction propensities. 
Simple polynomial rate laws can be provided as inline propensities and can 
be specified on the gateway level (e.g. in Matlab). 

The other option is to specify the rate laws in a model file written in C 
and recompile the Mex-interace. To simplify this process, we have supplied 
a makefile in '/msrc/make.m' using the mex command in Matlab with make 
options for two platforms (Linux and OS X). We have also provided an 
example GNU Makefile in '/src/Makefile' that allows full control of the make 
process. The example shows how to compile URDME on an Intel based 
MacBook using gcc and the Intel compiler ice. Both the Matlab and GNU 
Makefile will require some local editing, see Section 5.1. 

For other configurations, it should be fairly straightforward to modify 
the Makefiles in order to compile the code. When using Mex the user needs 
to make sure that the configuration is set up correctly and that a proper 
compiler is installed. In Section 5 we show examples of using both the inline 
propensities and compiled propensities. 

Note that one can easily use both inline and compiled propensities simul- 
taneously. This is convenient when only a few propensities are complicated 
and has to be compiled. 

4.2.1 Inline propensities 

An "inline propensity" is a compact data format for specifying basic chemical 
reactions with polynomial rate laws. An inline propensity P can be defined 
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Name 


Type 


Description 


1 V u X X o 


oCcJjicxi I illL J 


1 M Ulil U t/1 Ui o U U V (Jl Llilico . 


1 i o U c U X c; o 


oOdidi I llil J 


i>l UlllUcl Ui LlllltJitJIib b^JcClCb. 






X lilO diOU IXdiliCO lMIXUXO> — 






1 1 o ^ w -i_ ^ o /X XV \j ^ ^ ^ 1^ • 


Mypa ct i on ^ 


'^cpili^v ( 1 nt 1 


INTnTTinpr or rpapfiorm 


M1 
nx 


bCdldi I ill I J 


i>l UlllUcl Ui liiiiiit; JJl UjJt;iioltlc&. 


dsize 


scalar (int) 


Size of the data vector used in the 






propensity function. 


iiO 




LlVl / 1 gjlVOO UilC lillUlCli llLliliUd UI 




IllOL/Cv^XCO/SlMI.'CXXO / 


Gn PPl PG 7 in Gil R^TT^l 1 1 TTl P 1 
iDJJCV^lCO fc 111 O U. L/ V VJl LliliC 7* 








"CSPcLH 


Vt?L/LUl tU.ULlUlcJ 


xiii liiL/i cdolii^ ocL^UciiL/C Ui ^Uiiitb lii 






+1TV1P TxrViPfP trip G+Eitp riT tViP G^rGi'PTn ig 

tlillC WllClC tile OtdtC Ul tile oVOLClll lO 






tu uc i t; t Ui iicu. . 




Vector [Mreact ions-Ml] 


Propensity function pointers. See 




\I^X ^JLfCSllOX y ±r ULll / 


S^ppfinn 4 ^ for ^]p'^Pl^l'^i 

kJCl^UHJll '±.^ lUl LiCljClllO. 


report 


ReportFun 


Pointer to a report function. This 






function is called every time the 






chain reaches a value in tspan. 


vol 


Vector[Ncells] (double) 


The volume of the macroelements, 






i.e. the diagonal elements of the 






lumped mass-matrix M. 


sd 


Vector [Ncells] (int) 


The sub domain numbers of all sub- 






volumes. See Section 5.2 for more 






details. 


data 


Matrix (dsize xNcells) 


Generalized data vector. A pointer 




(double) 


to column j is passed as an additional 






argument to the propensities in sub- 






volume j. 



Table 4.1: The first few input arguments to rdme_solve. See also Table 4.2. 



9 



Name 



Type 



Description 



Sparse matrix (Ndofs) 
Ndofs) (double) 



Sparse matrix (Mspeciesx 
Mreactions) (int) 



Sparse matrix (Mreactions x 

[Mspecies+Mreactions]) 

(int) 



Matrix (3 x Ml) (double) 

Matrix (3 x Ml) (int) 

Sparse matrix (M^x Ml) (int) 
where Mg G [0, Ncells] 



The transpose of the diffusion matrix 
M-^K obtained from the FEM dis- 
cretization of the macroscopic diffu- 
sion equation, cf. (2.5). Each column 
in D (i.e. each row in M^^K) corre- 
sponds to a subvolume, and the non- 
zero coefficients D{i,j) give the dif- 
fusion rate constant from subvolume 
i to subvolume j. 

The stoichiometry matrix. Each col- 
umn corresponds to a reaction, and 
execution of reaction j amounts to 
subtracting the _7th column from the 
state vector. 

Dependency graph. The first 
Mspecies columns correspond to 
diffusion events and the following 
Mreactions columns to reactions. A 
non-zeros entry in element i of col- 
umn j indicates that propensity i 
needs to be recalculated if the event 
j occurs. See Section 5 for examples. 

The rate constants of the inline 
propensities. 

Specifies which species are involved 
in each inline propensity. 
Column j contains a list of all subdo- 
mains in which the jth inline propen- 
sity is off. For more details concern- 
ing the inline propensities, see Sec- 
tion 4.2. 



Table 4.2: Input arguments to rdme_solve (continued from Table 4.1). For 
more details, see the source file 'rdme.c'. For sparse matrices, the compressed 
column sparse (CCS) format is used. This is the same format Matlab uses 
and online documentation is available. 
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as 



^Ip + k2Xk + hn if i ^ j, 

^^^^ + k2Xk + hn iii = j. 

Here x is the column in x which contains the state of the subvolume con- 
sidered and Q is the corresponding volume. The coefficients and indices are 
specified in matrices K and I where K(:,r) = [ki, k2\ fca] and l(:,r) = [i\ j; k] 
are the constants corresponding to the rth inline propensity. The matrix S 
is a (possibly empty) sparse matrix such that S(:,r) lists all subdomains in 
which the rth inline propensity is turned off. Note that no inline propensities 
are active in subdomain zero! 




4.2.2 Compiled propensities 

The rate functions can also be supplied as an Ansi C model file. This is 
necessary if the rate functions are not simple polynomials or if they depend 
on additional input. Changing models in this setting amounts to recompiling 
the Mex gateway with the desired model file. 

Any model file must implement the following routines: 

• PropensityFun *ALLOC_propensities(void) 

• void FREE_propens it ies (PropensityFun *ptr) 

The first function should allocate and initialize an array of function pointers 
to the propensity functions and return a pointer to this array. The second 
function should deallocate the pointer ptr but sometimes additional actions 
need to be implemented. 

An important point is that neither Ansi C's malloc/free nor Matlab's 
mxMalloc/mxFree should be invoked explicitly. Instead it is preferred to use 
the defines MALLOC/FREE which can be modified at the stage of compilation 
(see for example '/msrc/make.m') 

The datatype PropensityFun is defined in the header 'propensities. h' 
(found in the 'include' directory) as 

typedef double (*PropensityFun) (const int *, double, 

const double *, int); 

Below is a commented example of a model file defining a simple chemical 
system composed of a single species X undergoing a dimerization reaction. 
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/* Propensity definition of a simple dimerization reaction. */ 
#include <stdlib.h> 
#include <stdio.h> 

/* Type definition of propensity functions: */ 
#include "propensities .h" 

/* Rate constant (in units M"{-I}s''{-1}) . */ 
const double k = le-3; 

double rFunKconst int *x, double vol, 

const double *data, int sd) 
/* Propensity for the reaction X + X -> 0. */ 
{ 

return k*x [0] * (x [0] -l)/vol; 

} 

PropensityFun *ALLOC_propensities (void) 

/* Allocation. */ 

{ 

PropensityFun *ptr = MALLOC(sizeof (PropensityFun) ) ; 
ptr[0] = rFunl; 

return ptr; 

} 

void FREE_propens it ies (PropensityFun *ptr) 

/* Deallocation. */ 

{ 

FREE (ptr) ; 

} 

More advanced examples can be found in the model files found in the 
'examples' directory; '/bistab/bistab.c' and '/Min/Min.c' corresponding to 
the examples discussed in Section 5. There we also explain how to compile 
and link the propensity functions. 

5 Examples 

In this section we provide a step-by-step description of how to set up and 
simulate two different models from the systems biology hterature. The first 
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is a model of a bistable system displaying phase separation when the species 
are diffusing slowly. This model was first studied in [5] for simple geometries. 
In this example, we will simulate the system in a model of a S. cerevisiae cell 
containing a nucleus and a large vacuole and we will see the ease with which 
we handle inner (curved) boundaries. 

The second example is a model of protein oscillations in a model of an 
E. coli cell proposed in [12]. Here we will illustrate how to use different 
subdomains in order to incorporate membrane diffusion in the model. 

5.1 Installing URDME 

There is no automated installation procedure. Simply download the source 
and decompress it in a directory of your choice (it will result in the path 
to the code being '/yourdir/urdme/'). There are a few path dependencies 
that need local editing. In the following, all paths are given relative to 
' / yourdir / urdme / ' . 

1. If you plan to use the GNU Makefile, open it ('/src/Makefile') and edit 
the variables PREFIX, MATLAB_EXT and MATLAB_LIB. In addition, you 
may need to modify other variables, depending on your configuration. 

2. Open the file '/msrc/make.m'. In Matlab, use the function mexext to 
find out the extension of Mex-files on your platform. If it is something 
different from the provided examples ('mexa64' and 'mexmaci'), you 
will have to add an option similar to the existing ones. The compile 
and link flags of the first example ('mexa64') are likely to work for most 
Linux fiavors and for Solaris (using gcc). The second option will likely 
work for a wide range of OS X configurations. More information on 
how to build a valid Mex-file on your platform can be obtained from 
the file 'mexopts.sh' (or 'mexopts.bat' in a Windows environment) in 
your Matlab installation directory. 

Before proceeding, make sure that URDME compiles. In Matlab, change the 

directory to '/msrc' and type 

>> make . ./examples/bistab/bistab.c 

or using the GNU Makefile, from a bash-shell, change the directory to '/src' 
and type 

$ make MODEL=. ./examples/bistab/bistab 

If mexrdme was built successfully, you can proceed to the examples in the 
two following sections. 

Currently, there is no support for Windows. To use URDME in a Win- 
dows environment we recommend installing gcc for Windows (such as in 
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Mingw) and using for example the freely available tool gnumex to configure 
mex to use gcc. 

5.2 General workflow 

In this section we describe the general workflow involved in setting up and 
simulating a model using the Comsol and Matlab interfaces. In the following 
sections, we then consider two specific examples in some detail. 

1. Specify the model. In order to do this, you need to define the ge- 
ometry as well as the chemical reactions in the model. In Comsol 
Multiphysics 3.4 a geometry can be created and a discretization of 
the diffusion equation with Neumann boundary conditions is readily 
obtained. In Matlab, a model file that defines the stoichiometry ma- 
trix N and the dependency graph G can easily be specified. Depending 
on how the reactions are to be given, this file may also define the inline 
propensities. 

2. Convert the FEM structure to valid input. After meshing your model, 
export the FEM structure to Comsol Script and create a model .mat 
file for use by the Matlab interface. This can be achieved by calling 
the interface routine f em2rdme. 

3. Run the simulation. Load the model file in Matlab, specify initial condi- 
tions and (if you are using compiled propensities) compile the gateway 
mexrdme using either 'make.m' or the GNU makefile, cf. Section 5.1. 
Then run the simulation by calling rdme_solve via the interface rdme. 

4. Postprocessing. Save the result to a result .mat file. In Comsol Script, 

load this file and add the solution to the FEM struct previoTisly ex- 
ported by calling rdme2f em. At this point, you can use all postpro- 
cessing options available in Comsol. If you have other needs not cov- 
ered by the built-in routines, you can implement your own routines in 
e.g. Matlab. 

5.3 A bistable system in S. cerevisiae 

In this section we will simulate the model from [5] in a spherical volume 
modeling a Brewer's yeast cell. It contains a nucleus and a large vacuole, 
which will be treated as solid objects into which no molecules can enter. The 
set of chemical reactions in the system can be found in Table 5.1. All the 
files needed to complete this example can be found in the directory '/exam- 
ples/bistab/'. 
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Ea^Ea + A Eb^Eb + B 

Ea + B^ EaB Eb + A^ EbA 

kd kd 

EaB + B^EaB2 EbA + A^EbA2 
kd kd 

Table 5.1: The chemical reactions of the bistable model. The constants take 
the values h = 150s-^ ka^ 1.2 x lO^s^^M-^ ka = lOs^^ and k4 = 65"^ 



1. Start Comsol Multiphysics 3.4. In the Model navigator, select 'Open' 
and open the file 'bistab.mph'. This model file specifics the yeast ge- 
ometry and in the Draw mode you can study the different subdomains. 

2. The model is composed of three subdomains. We will treat the inner 
structures as solid objects that molecules cannot enter. This can be 
achieved in different ways. One of them is to set the diffusion coefficient 
identically to zero in these subdomains. Here we will remove these inner 
structures and create holes in the geometry which reduces the number 
of degrees of freedom. In the Draw mode, select all three subdomains 
by shift-clicking on them in turn. In the 'Draw' menu, select 'Create 
Composite Object' and choose 'Difference' (this can also be done using 
the action buttons in the draw mode action bar). Next, again in the 
'Draw' menu, select 'Delete interior boundaries'. 

3. Next we will initialize the mesh. In the Mesh menu, choose 'Free Mesh 
Parameters'. Set the maximum element size to 0.1 x 10~^m and click 
on the 'Remesh' button. 

4. The parameters in the model (in this case the diffusion coefficients of 
the different species) can be edited in the menu 'Physics - Subdomain 
Settings'. Choose subdomain one and check that the species all have 
the diffusion constant 7 = 0.5 x 10~^^m^/s. 

5. We are now ready to export the FEM- structure to Comsol Script in 
order to prepare for the simulation. In the 'File' menu, choose 'Export 
- FEM Structure...' and choose the name of the struct, e.g. 'fem'. This 
opens the Comsol Script window and you can display the contents in 
'fem' by simply typing 

C» fem 

6. We have now specified the geometry and will next create the input to 
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rdme_solve. First, change the working directory to '/msrc' and set the 
paths by caUing startup. 

C» cd /yourdir/urdme/msrc ; 
C» startup; 

In order to assemble the input to rdme_solve, call the interface routine 
fem2rdme, 

C» fem2rdme(fem, 'bistab') ; 

This will generate the file 'bistab.mat' that contains the diffusion ma- 
trix, the subvolume sizes as well as information regarding the different 
subdomains (in this example there is only one). 

7. Start Matlab, change the working directory to '/msrc' and set the paths 
as above. Then load the model file we just created 

» load bistab.mat; 

and examine the variables. We will now specify the set of chemical re- 
actions. There are two different ways of doing this; using inline propen- 
sities or using compiled propensities, see Section 4. Here we will use 
the compiled propensities specified in the model file 'bistab.c'. Open 
the file and look at the reaction definitions. Compile the code with 
these propensities 
» make . ./examples/bistab.c; 

8. Before we can run the simulation, we need to specify the stoichiomctry 
matrix N and the dependency graph G. They are both defined in the 
script 'bistab.m'. Run it by invoking 

» bistab; 

9. We have yet to specify the initial conditions. Let the enzymes Ea and 
Eb be present in a total of 100 molecules, randomly distributed in the 
whole domain. The other species are set to zero everywhere. Something 
like this can be achieved by 

>> indl = randperm(Ncells) ; 

>> ind2 = randperm(Ncells) ; 

» u0(3,indl(l:100))=l; 

» u0(4,ind2(l:100))=l; 

10. At this point, we are ready to run the simulation by invoking rdme_solve 
through rdme; 

» bistab.sol = rdme(0:5: 1000,uO,D,N,G,vol,sd, [2 123]); 
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When the simulation has finished, save the solution and the volume 
vector (it is needed for the postprocessing in order to scale the solution 
to local concentrations). 

» save('bistab_results' , 'bistab_sol ' , 'vol'); 

11. In Comsol Script, load the solution 
C» load bistab_results .mat 

and add the solution to the FEM-struct using the interface routine 
rdme2f em; 

C» fem = rdme2fem(bistab_sol, vol, 0:5:1000); 

12. In the main Comsol window, import the FEM structure that now con- 
tains the solution. This is done in the File menu; 'File - Import - FEM 
Structure...'. Change mode to 'Plot mode' and look at the solution. 
You can change variables and the type of plot in the menu 'Postpro- 
cessing - Plot Parameters...', where you also have the possibility to 
make animations. This can also be done directly in Comsol Script 
using the functions postplot and postmovie, see the documentation 
for more details. 

Figure 5.3 shows an example of the output from a simulation of the model. 
For more details concerning the interpretation of the results, see [5, 10]. 

5.4 MinD/MinE oscillations in E. coli 

In this section we will reproduce simulations of the MinD /MinE system stud- 
ied in [12] in a model of an E. coli bacterium. It is rod shaped with length 
3.5|xm and diameter 0.5|j.m. The reactions and parameters of the model can 
be found in Table 5.2. 

MinD_c_atp MinD_m MinD_c_atp + MinD_m 2MinD_m 
Min_e+MinD_m ^ MinDE MinDE MinD_c_adp + Min_e 

k 

MinD_c_adp MinD_c_atp 

Table 5.2: The chemical reactions of the MinD/MinE model. The constants 
take the values kd = 0.0125/im~^s~\ kdo = 9 x 10^M~^s"^, kde = 5.56 x 
WM-^s-\ ke = Q.7s-^ and kp = 0.5s-^ 

Two of the species, MinD_m and MinDE are bound to the membrane and 
can not diffuse in the cytoplasm. In order to incorporate this in the model 
we need to define two different subdomains, the membrane and the cytosol. 
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V . 



(a) t = 250s 




(c) t = 900s 



(b) t = 250s 



■ ' 1 



(d) t 900s 



Figure 5.1: Snapshots of the system simulated with 7 = 0.5 x IQ^^^m? / s. 
The left column shows species A and the right species B. They are present 
in different parts of the geometries, and the global behavior does not display 
bist ability. 
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and specify which reactions take place in each subdomain. We also need 
to modify the diffusion matrix in order to inactivate diffusion events from 
the membrane into the interior of the cell. This is achieved by using the 
subdomain data vector sd. 

Note that the subdomain data vector sd is a completely general piece of 
information. In particular, it need not he explicitly related to the 'subdomain' 
numbers in Comsol. 

One of the reactions taking place on the membrane has a scaling with a 
local length scale. In this example, we will use inline propensities to which 
we pass the average of the side length of the tetrahedra on the boundary. 
If compiled propensities were used, it would instead be possible to pass the 
actual length of each individual subvolume to the propensity functions using 
the generalized data matrix data. All the files used in this example can be 
found in the directory ' /examples/ Min'. 

1. Open Comsol and select 'Chemical engineering module - Mass trans- 
port - Diffusion - Transient analysis' (3D). In the 'Dependent variables' 
field write MinD_c_atp, MinD_m, Min_e, MinDE, MinD_c_adp. These 
are the names of the variables that we will use. Select Lagrange - 
Linear elements and press 'OK'. 

2. Next we create the geometry. Wc will build the rod shaped domain 
from two spheres and one cylinder. Press the 'Cylinder' button and in 
the radius and height field enter 0.5e-6 and 3.5e-6 and press 'OK'. 
You should now see a cylinder in your workspace. In the 'Draw mode' 
action bar, press 'Zoom extents' in order to get a better view of the 
domain. Press the 'sphere' button and enter . 5e-6 in the radius field 
and press 'OK'. Create another identical sphere but enter 3.5e-6 as 
the z-coordinate. Select all three figures and press the 'union' button 
and then the 'Delete interior boundaries' button. 

3. Having defined the geometry, the next step is to specify the parameters 
in the model. In the menu 'Physics - Subdomain settings', choose sub- 
domain 1 and set the diffusion constants to 2.5e-12 for MinD_c_adp, 
MinD_c_atp and Min_e. For MinDE and MinD_m the diffusion constant 
should be le-14. 

4. We must also create two domains. One interior domain that represent 
the cytoplasm and one boundary domain that represent the membrane. 
This is done by defining rdme_sd as an expression with different value 
in the different subdomains. It can then be used to find the nodes on 
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the boundary and in the interior. Select 'Options - Subdomain expres- 
sions' and enter rdme_sd with value 1 and click 'OK'. Select 'Options 
- Boundary expressions' and select all boundaries (there should be 12 
of them). Enter rdme_sd with value 2. Finally select 'Options - Global 
expressions' and enter rdme_sdlevel with value 2 indicating that the 
lowest dimension where rdme_sd is defined was on the surfaces. 

5. In the 'Mesh' menu, select 'Mesh - Free mesh parameters' and choose 
'Custom mesh size'. Set the maximum element size to . 7e-7 and press 
'Initiahze mesh'. Now select 'Solve - Update model' and export the 
FEM structure from the 'File - export' menu. This will open Comsol 
Script and the FEM structure will be available in the workspace. 

6. The information we need regarding the discretization of the geometry 
is stored in the FEM structure (call it f em). At this point, we need to 
transform the information in fern into the input format of rdme_solve. 
This is done by calling the interface routine f em2rdme. First, initialize 
URDME by caUing startup (from the '/msrc' directory). Then, as- 
semble the input and store it in the file 'min.mat' by 

C» f em2rdme(f em, 'min' , 1) ; This will save the variables uO, D, vol, 
data and sd to the file min.mat. In order to make sure that the diffu- 
sion of the species MinD_m and MinDE are active only on the boundary, 
type the following: 

C» load min; 

C» icells = f ind(sd==l) ; 

C» D = min_membrane(D, icells, [2 4], 5); 

C» save ('min' , 'tspan' , 'uO' , 'D' , 'vol' , 'data' , 'sd' , '-mat') ; 

This will set the diffusion rate constants to zero in the interior. The 
first command loads the variables that fem2rdme saved. The second 
command gives the indices to the cells in subdomain one, which we 
defined to be the interior of the domain by the expression rdme_sd. The 
third one calls the script minjnembrane which turns off the diffusion 
for species 2 and 4 (MinD_m and MinDE). Finally, we rewrite the file 
min.mat with the modified matrix D. 

7. We are now ready to specify the part of the model related to the chemi- 
cal reactions. Open Matlab and initialize URDME by calling startup. 
Load the file min.mat. Next, in a new Matlab-script create the stoi- 
chiometric matrix and the dependency graph which are used to update 
the state due to chemical reactions and to minimize the calculation of 
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reaction rates respectively. If you want to skip this step, just use the 
script Min.m, where they are aheady defined. 



N = 



1 

-1 







1 

-1 









1 

1 -1 

-1 1 

-1 



-1 







1 



G 



1 1 1 1 

1 1 1 1 1 1 

110 11110 

1 1 1 

1 1 1 



We will use inline propensities so we also need to create K, I, and S 



K 








kdD 

















kp 




Note the unit of kd, it involves the local length scale of the computa- 
tional cells at the boundary. Using inline propensities, we cannot pass 
this to the propensity functions. Instead, we use the average length 
scale and modify accordingly. The local length scale can be ob- 
tained using the Comsol routine posteval and for this mesh size the 
average at the boundary takes the value 0.041 x 10~^. 



113 11 
12 2 11 
1114 5 

S = [ 1 1 1 1 ] . 

We have now specified both the geometry, the chemical reactions and 
the dependency graphs. Before we run the simulation, we need to 
set the initial condition. Set the initial number of molecules to 2001 
MinD_c_atp, 2001 MinD_c_adp and 1040 MinE, randomly distributed. 



» iiiind=[ceil(rand(2001, l)*Ncells)] ; 
» mind2=[ceil(rand(2001,l)*Ncells)] ; 
» mine=[ceil(rand(1040, l)*Ncells)] ; 
» a=spcirse(mind, 1 , 1 ,Ncells , 1) ; 
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» b=sparse(mine, 1, l.Ncells, 1) ; 

» c=sparse(mind2,l,l,Ncells,l) ; 

» uO=[a zeros (Ncells.l) b zeros (Ncells, 1) c] ' ; 

» uO=full(uO) ; 

9. Set tspan= [0 : 30] and start the simulation by: 

» U = rdme (tspan,uO,D,N,G, vol, data, sd, [2 123], K, I, S) ; 

After the simulation has finished, save the solution on a .mat file, e.g. 
» save('min_results' , 'U') ; 

10. To visualize the trajectory, load the result in the same Comsol Script 
window used before. Make sure that the variable vol is in your workspace 
and convert the result to Comsol's internal representation using the in- 
terface routine rdme2f em; 

» fem = rdme2fem(fem,U,vol, [0:30]) ; 

11. Import the FEM structure into the Comsol GUI from the 'File - Import' 
menu and plot the solution in the postprocessing mode. For example of 
output generated by this model, see Figure 5.2, where we have plotted 
the variable MinD_m on the membrane using 'Boundary Plot' from the 
'Plot Parameters' menu. 



6 Performance 

The purpose of this section is to give an indication of the performance of the 
code by simulating the model in Section 5.3 on increasingly fine meshes. The 
time required to simulate a particular model will obviously depend on the 
parameters since the average timestep is the inverse of the sum of all reaction- 
and diffusion rates. For example, the fraction of diffusion to reaction events 
increases with the resolution and the number of events increase linearly with 
the diffusion constant. Consequently, it is hard to tell in advance how many 
events that are required to reach a given final time. In Table 6.1 we show for 
fixed parameters in the model how the number of events simulated per hour 
depends on the number of subvolumes. 

The geometry was taken to be a cube with side length 6 x 10~^m. At 
t — Os, there were 300 molecules each of Ea and Eb, randomly distributed 
and all other species were set to zero. The final time was chosen to be 
T = 300s, and the state was saved every three seconds. The timings in 
Table 6.1 does not include the time taken to generate the mesh and the 
input to URDME, which was much less than the simulation time. 
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(a) t = 10s 



(b) t = 40s 




(e) t ^ 65s (f) t = 70s 

Figure 5.2: Simulation of the MinD model in E. coli. At t = Os, no MinD 
is present on the membrane. After an initial period when the number of 
molecules are increasing over the whole membrane in (a)-(b), the membrane 
bound MinD oscillates from pole to pole in (c)-(f). 
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Ncells 


11,107 


25,608 


83,843 


CPU time [s] 


155 


230 


667 


% diffusion events 


55.7 


69.1 


83.3 


events/h [xlO^] 


2.04 


1.95 


1.24 



Table 6.1: Performance when the spatial resolution is increased. Note that 
the fraction of diffusion events increases with the resolution of the mesh. 

The simulations were made on a Macbook Core Duo 2.0 Ghz with 2GB 
RAM, using compiled propensities and compiled using gcc 4.3.2 with op- 
timization flags -03 -ftree-vectorize. 
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A Algorithms 



One of the most popular algorithms to generate realizations of the CTMC in 
the well-stirred case is Gillespie's direct method (DM) [15]. Several algorith- 
mic improvements of this method exist, one of them being the next reaction 
method (NRM) due to Gibson and Bruck [14]. 

The underlying algorithm in URDME is the next subvolume method 
(NSM) [5]. The NSM essentially combines ideas from NRM and DM in 
order to tailor the algorithm to reaction-diffusion processes. 

For reference, we first state below both DM and NRM and then outline 
NSM. 

Algorithm 1 Gillespie's direct method (DM) 

Initialize: Set the initial state x and compute all propensities a;r(x),r = 
1, ■ ■ ■ , ^reactions- Also set t = 0. 
while t<T do 

Compute the sum A of all the propensities. 

Sample the next reaction time (by inversion), r = — log(rand)/A. Here 
and in what follows, 'rand' conveniently denotes a uniformly distributed 
random number in (0, 1) which is different for each occurrence. 
Sample the next reaction event (by inversion); find n such that 

Ej^i^^iW < A rand < 

Update the state vector, x = x — iV„ and set t = t + t. 
end while 
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Algorithm 2 Gibson and Bruck's next reaction method (NRM) 

Initialize: Set t — and assign the initial number of molecules. Generate 

the dependency graph G. Compute the propensities cj,.(x) and generate 
the corresponding absolute waiting times for all reactions r. Store those 
values in a heap H. 
while t<T do 

Remove the smallest time t„ = Hq from the top of H, execute the nth 
reaction x := x — and set t Tn- 
for all edges n — > j in G do 
if J 7^ n then 

Recompute the propensity uij and update the corresponding waiting 
time according to 

^old 

new^^^/ old_ 

3 \ 3 I , .new 

else {j = n\ 

Recompute the propensity a;„ and generate a new absolute time 
. Adjust the contents of B. by replacing the old value of t„ with 
the new one. 
end if 
end for 
end while 
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Algorithm 3 The next subvolume method (NSM) 

Initialize: Compute the sum (t[ of aU reaction rates Uri and the sum af of 
all diffusion rates OjjX^j in all subvolumes i = 1, . . . , N^ells- Compute the 
time until the next event in each subvolume, Xj = — log (rand)/ ((t[ + af), 
and store all times in a heap H. 
while t < T do 

Select the next subvolume (n where an event takes place by extracting 
the minimum r„ from the top of H. 
Set t = Tn- 

Determine if the event in C,n is a reaction or a diffusion event. Let it be 
a reaction if (crj^ + cr^) rand < crj^, otherwise it is a diffusion event, 
if Reaction event then 

Determine the reaction channel that fires. This is done by inversion 

of the distribution for the next reaction given r„ in the same manner 

as in Gillespie's direct method in Algorithm 1. 

Update the state matrix using the (sparse) stoichiometry matrix N. 

Update cr^j and 0"^^ using the dependency graph G to recalculate only 

affected reaction- and diffusion rates, 
else {Diffusion event} 

Determine which species Sin diffuses and subsequently, determine to 

which neighboring subvolume Cn'- This is again done by inversion 

using a linear search in the corresponding column of D. 

Update the state: Sni = Sni-l, Sn'i = Sn'i + 1. 

Update the reaction- and diffusion rates of subvolumes Cn and Cn' using 

G. 
end if 

Compute a new waiting time t„ by drawing a new random number and 
add it to the heap H. 
end while 
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